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SUMMARY 

A  comparison  between  the  experimental  visualization  and  numerical 
simulations  of  the  occurrence  of  vortex  breakdown  in  laminar  swirling  flows 
produced  by  a  rotating  endwall  is  presented.  The  experimental  visualizations 
of  Escudier  (Experiments  in  Fluids,  2,  1984)  were  the  first  to  detect  the 
presence  of  multiple  recirculation  zones  and  the  numerical  model  presented 
here,  consisting  of  a  numerical  solution  of  the  unsteady  axisym metric  Navier- 
Stokes  equations,  faithfully  reproduces  these  newly  observed  phenomena  and 
all  other  observed  characteristics  of  the  flow.  Part  II  of  the  paper  examines 
the  underlying  physics  of  these  vortex  flows. 
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1.  Introduction 

Vortex  breakdown  in  swirling  flows  has  been  the  subject  of  much  attention  since  it 
was  first  recognized  in  the  tip  vortices  of  delta  winged  aircraft  (Peckham  &  Atkinson 
1957).  More  recently  it  has  emerged  as  a  serious  problem  at  high  angle  of  attack  for 
highly  manoeuvrable  military  aircraft.  As  yet,  there  is  no  general  consensus  on  the  phys¬ 
ical  mechanisms  responsible  for  its  occurrence.  The  term  vortex  breakdown  is  associated 
with  an  abrupt  change  in  the  character  of  a  columnar  vortex  at  some  axial  station.  It 
is  usually  observed  as  a  sudden  widening  of  the  vortex  core  together  with  a  deceleration 
of  the  axial  flow  and  is  often  followed  by  a  region  or  regions  of  recirculation. 

Experimental  studies  of  vortex  breakdown  over  delta  winged  aircraft  have  been 
severely  hampered  by  the  sensitivity  of  the  flow  to  the  presence  of  external  probes, 
making  quantitative  measurements  of  the  flow  in  the  ‘burst’  regions  much  more  difficult. 
Further  complications  arise  as  a  result  of  the  large  number  of  parameters  involved, 
some  of  which  are  difficult  to  measure  and/or  control  and  whose  importance  to  the 
mechanisms  leading  to  vortex  breakdown  is  not  known.  With  the  identification  of  tlm 
vortex  breakdown  of  swirling  flows  in  cylindrical  tubes  (Harvey  1960),  a  large  number  of 
experimental  investigations  were  undertaken  (eg.  Sarpkaya  1971;  Faler&  Leibovich  1978; 
Escudier  Sz  Zehnder  1982).  A  number  of  distinctive  forms  of  breakdown  were  observed 
in  these  investigations,  which  also  showed  the  phenomenon  typically  to  be  unsteady  and 
possessing  various  degrees  of  turbulence  and  asymmetry.  From  these  results,  it  had  been 
assumed  that  a  description  of  vortex  breakdown  would  require  a  three-dimensional,  time 
dependent  calculation  of  the  Navier-Stokes  equations  albeit  with  the  possibility  at  high 
Reynolds  numbers  of  modelling  the  small  scale  turbulence. 

Recently,  Escudier  (1984)  observed  the  phenomenon  of  vortex  breakdown  in  swirling 
flows  in  a  cylindrical  container  with  a  rotating  endwall,  using  a  laser-induced  fluores¬ 
cence  technique.  The  experimental  results  extended  those  obtained  earlier  by  Vogel 
(1968)  and  Ronnenberg  (1977)  and  are  the  first  tc  be  reported  in  which  multiple  break¬ 
down  bubbles  exist  in  the  closed  cylindric.  :  coo-  ,<  try.  The  recirculation  bubbles  were 
observed  to  be  axisymmetric  and  steady  over  irge  range  of  the  governing  parame¬ 
ters,  which  re  inforced  Escudier’s  (1981)  view,  earlier  expressed  as  a  result  of  a  series 
of  swirling  pipe  flow  experiments  (Escudier  &•  Keller  1983),  that  vortex  breakdown  in 
general  is  inherently  axisymmetric  and  that  departures  from  axial  symmetry  result  from 
instabilities  n ■  >t  directly  associated  with  the  breakdown  process.  Escudier’s  experiments 
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provide  a  particularly  well-defined  set  of  flows  in  which  vortex  breakdown  occurs  and 
his  figures  (1984)  motivated  the  present  study  in  which  the  Navier-Stokes  equations  are 
solved  numerically  under  the  assumption  of  axial  symmetry.  A  principal  aim  of  this 
study  was  to  clarify  whether,  under  the  assumption  of  axial  symmetry,  the  numerical 
solutions  accurately  reproduce  all  of  the  observed  flow  characteristics  reported  by  Es- 
cudier  (1984).  A  further  motivation  for  undertaking  the  investigation  was  to  develop  a 
more  detailed  understanding  of  the  physics  of  the  flow  and  to  clarify  features  that  were 
not  readily  resolved  from  the  visualizations.  For  example,  the  dye  streaks  are  ‘blurred’ 
in  the  region  of  critical  points  in  the  experiment  and  the  nature  of  the  downstream  end 
of  the  recirculation  bubble  is  difficult  to  ascertain  from  flow  visualization.  A  main  reason 
for  this  ‘blurred’  picture  is  dye  diffusion  but,  in  some  cases  Escudier  (1984)  noted  that 
certain  recirculation  zones  never  develop  a  well  defined  internal  structure.  The  numeri¬ 
cal  study  shows  that  these  regions  in  fact  possess  a  localized  low  frequency  unsteadiness 
and  this  unsteadiness  would  result  in  a  'blurred'  flow  visualization. 

In  part  II  of  this  study  the  physics  of  the  flows  is  examined  in  more  detail.  The 
solutions  are  considered  from  the  point  of  view  of  standing  centrifugal  waves  and  an 
understanding  of  the  mechanism  responsible  for  the  vortex  breakdown  of  these  con¬ 
fined  swirling  flows  is  presented.  The  mechanics  of  vortex  breakdown  in  this  particular 
geometry  suggests  a  generalization  which  is  discussed  for  some  other  geometries. 

Previous  numerical  studies  of  swirling  flows  confined  within  a  cylindrical  container 
(eg.  Pao  1970;  Lugt  &  Haussling  1973;  Dijkstra  &  Heijst  1983)  have  either  been  restricted 
to  small  aspect  ratio  flows  ( H/R ,  ratio  of  height  and  radius  of  the  cylinder,  less  than 
1),  or  flows  with  Re,  the  rotational  Reynolds  number,  too  small  for  vortex  breakdown 
to  occur.  Lugt  &  Haussling  (1982)  is  a  notable  exception  in  which  H/R  <  16.  and  the 
occurrence  of  a  single  breakdown  bubble  is  reported.  Their  computational  mesh  does 
not  appear  to  be  fine  enough,  however,  since  they  found  the  size  and  location  of  the 
bubble  to  vary  as  the  resolution  was  altered  and  the  time  scales  of  their  solutions  are 
not  consistent  with  experimental  observations. 

Since  submitting  the  first  version  of  this  paper,  Lugt  and  Abboud  (1987)  have 
also  published  numerical  computations  of  the  vortex  breakdown  phenomenon  produced 
in  an  enclosed  cylinder  by  a  rotating  cndwall.  Their  calculations  and  those  presented 
here  overlap  somewhat  as  both  studies  use  the  available  experimental  results  of  Escudier 
(1984)  as  comparisons  for  the  numerical  calculations.  The  results  presented  by  Lugt  and 
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Abboud  (1987)  generally  agree  very  well  with  those  presented  here  and  their  numerical 
method  of  integrating  the  governing  equations  is  quite  different  to  the  present  method. 
A  large  number  of  their  results  had  not  been  integrated  forward  until  a  true  steady 
state  had  been  reached  and  no  information  regarding  the  distribution  of  the  azimuthal 
component  of  vorticity  is  presented,  which  will  be  shown  in  Part  II  of  this  study  to  be 
crucial  in  understanding  the  physical  mechanism  responsible  for  the  vortex  breakdown 
of  the  vortical  core  flow. 


4 


2.  Governing  Equations  and  Boundary  Conditions 

The  flows  of  interest  are  simulated  by  considering  a  cylinder  of  radius  R  and  length 
H,  the  bottom  endwall  of  which  is  impulsively  started  to  rotate  at  constant  angular  ve¬ 
locity  Q.  The  fluid,  which  completely  fills  the  cylinder,  is  incompressible,  of  uniform  den¬ 
sity  and  a  constant  kinematic  viscosity  v.  The  axisymmetric  form  of  the  Navier-Stokes 
equations,  in  cylindrical  co-ordinates  (r,<p,z)  with  corresponding  velocity  components 
(u,u,u>),  is  employed.  Time  and  length  are  scaled  by  fl  and  l/R  respectively. 

The  system  of  equations  is  solved  by  employing  the  streamfunction-vorticity  formu¬ 
lation,  where  the  pressure  does  not  appear  explicitly.  This  is  achieved  by  the  introduction 
of  a  streamfunction  ip,  where 

dip  ,  1  dip  . 

u  —  a-  and  10=-  —  ,  (1) 

r  dz  r  or 

which  satisfies  continuity  and  gives  for  the  azimuthal  component  of  vorticity 

\d2ip  _  9  (i^t\ 

^  r  dz2  dr  r  dr 

Incorporating  (1)  and  (2)  in  the  Navier-Stokes  equations  leads  to  the  following  predic¬ 
tion  equations  for  the  azimuthal  components  of  velocity  and  vorticity  together  with  the 
prognostic  equation  for  the  streamfunction: 


dv  dv 

n  =  J(r)  +  2:5i  + 


1  .d2 v  d2v  1  dv  v  . 

Re  dz 2  dr 2  r  dr  r2 1  ’ 

(3) 

J,  cPj  dS  \  drj 

Re  dz 2  dr2  r  dr  r2  ’ 

(4) 

p  1  dip 

(5) 

2  rTr  ~  rVl 

dip  d  dip  d 
dz  dr  dr  dz 

The  boundary  conditions  relevant  to  the  experiments  of  Escudier  ( 1984)  to  complete 
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the  system  (3)-(5)  are: 

ip  =  v  =  7/  —  0  (r  =  0,  0  <  z  <  H / R), 

ijj  =  v  =  0,  7  = -"Try  (r  =  1,  0  <z<H/R), 

r  Or1 

1  (®) 

=  0,  v  —  r,  V  =  ---^~2  (z  ~  °>  0  <  t-  <  1), 

r  f72z 

jI>  =  v-0,  n  =  -~TTJ  (z  ~  H/R>  0<r<l). 
r  ozL 

The  boundary  condition  at  r  =  0  is  due  to  the  axial  symmetry  of  the  flow,  the 
boundaries  at  z  =  H/ R  and  r  =  1  are  rigid  and  stationary,  while  at  z  =  0,  the  rigid 
endwall  is  in  constant  rotation  for  t  >  0.  The  fluid  is  stationary  for  t  <  0. 
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3.  Method  of  Solution 

The  system  of  differential  equations  (3)  -  (5)  is  replaced  by  an  approximating  set 
of  finite-difference  equations,  defined  on  a  uniform  mesh  in  the  r,  z  plane  at  times  nAt 
(n  =  0, 1,2, ...).  There  are  (nr  +  1)  x  (nz  +  1)  grid  nodes,  and  in  all  the  cases  presented 
here,  the  increments  Ar  and  A z  are  equal.  Second-order  accurate  centered  differences 
are  used  to  approximate  all  the  spatial  derivatives  except  those  in  the  advection  terms, 
where  the  Jacobian  conserving  difference  operator  J  is  approximated  using  a  scheme  due 
to  Arakawa  (1966).  The  time-differencing  scheme  used  is  that  of  Miller  &  Pearce  (1974) 
which  consists  of  alternate  timesteps.  For  a  prediction  equation  of  the  form  dF/dt  - 
G(F),  the  odd  time-steps  are  advanced  according  to 

F’  =  Fn  +  A  tG(Fn), 

Fn+i  =  Fn  +  AlG(F-)' 

while  for  even  time-steps 

Fn+l  =  pn  +  A  lG(Fny 

The  prognostic  equation  (5)  for  the  streamfunction  is  solved  using  the  generalized  cyclic 
reduction  method  of  Sweet  (1974). 

To  implement  the  boundary  conditions  for  the  azimuthal  vorticity  at  the  solid 
boundaries  at  time  (n  +  l)At,  the  streamfunction  field  at  that  time  is  needed  This 
is  achieved  by  the  following  algorithm.  First,  the  prediction  equations  (3)  and  (4)  for  t' 
and  7/  may  be  advanced  to  the  next  time-step  (n  -f  1)A(  for  the  interior  points  since  the 
prognostic  equation  (o)  fur  ip  only  requires  knowledge  of  the  interior  values  of  r/.  From 
the  solution  to  (5)  together  with  the  boundary  conditions  for  ip,  the  streamfunction  is 
known  everywhere  at  t  =  (n  4-  l)At.  The  boundary  values  of  tj  can  now  be  estimated 
from  xp  by  noting  that  ip  and  its  normal  derivative  vanish  at  a  rigid  boundary.  Hence, 
expanding  ip  about  the  first  point  in  from  the  boundary  leads  to 

7](nr  +  \  ,j)  -  -2ip(nr,j)/ Ar2 , 

V(i,  1)  =  -2^(.',2)/(rAz2), 

and 

-1)  ~ 2V  (t,  nz)/(rAz2). 
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Proper  resolution  of  the  boundary  layers  and  any  recirculation  zones  was  ensured 
by  the  use  of  A r  —  Ac  =  1/60.  A  number  of  runs  for  H/R  =  1.5  were  performed 
using  Ar  =  A z  —  1/40  with  only  slight  differences  in  the  outline  of  the  recirculation 
zone  being  noted  —  the  positions  of  the  stagnation  points  were  not  appreciably  altered 
by  the  reduction  in  resolution.  One  run  at  H/R  —  2.5  and  Re  =  2126  implemented 
Ar  =  A z  =  1/100.  The  difference  between  the  streamlines  from  these  results  and  those 
from  calculations  with  Ar  =  A z  =  1/60  were  insignificant  All  the  calculations  presented 
here  are  for  Ar  =  A z  =  1/60  and  At  =  0.05,  except  for  the  one  case  of  H/R  =  2.5 
and  Re  =  2126  where  A r  =  Az  =  1/100  and  At  =  0.025.  In  all  cases,  the  stability 
of  the  model  was  ensured  by  the  time-step  satisfying  both  the  Courant-Friedrichs-I.ewy 
condition  and  the  diffusion  requirement  (Williams  1967), 

1  ^ 

At  <  -ReAr. 

f  The  quality  of  the  numerical  computations  was  further  verified  by  (i)  direct  com¬ 

parisons  with  experimental  results  (§4),  and  by  (ii)  evaluation  of  the  vorticity  integral 
(Dijkstra  &  van  Heijst  1983).  From  the  definition  of  the  azimuthal  vorticity  (2),  together 
with  the  boundary  conditions,  the  streamfunction  and  its  normal  derivatives  vanish  on 
solid  boundaries  and  4>  -  0  on  the  axis  of  symmetry,  hence 

//(Vl;  =  2'£o£"r!^r  =  0.  (7) 

At  t  =  1 000,  this  integral  was  evaluated  by  means  of  a  two-dimensional  trapezoidal 
rule,  and  the  positive  and  negative  contributions  to  the  trapezoidal  sum  accumulated 
separately  in  /(q^)  and  /(q  )  respectively.  In  the  limit  Ar,  Az  — >  0,  /(q  +  )  +  /(q  )  <>. 

The  degree  to  which  this  is  achieved  is  an  indication  of  the  quality  of  the  numerical 
scheme.  Table  1  lists  the  values  of  the  integral  (7)  together  with  /(q  +  )  and  / ( q '  )  for 
a  selection  of  the  cases  considered.  It  is  clear  that  the  accuracy  of  the  calculation* 
is  improved  by  increasing  the  resolution  from  Ar,  Az  =  1/40  to  1/60.  For  the  case 
Re  =  1256,  lTt i  =  (/(q~)  +  ~  /(q+))  —  0.0015  when  Ar,  Az  =  1/40, 

whereas  Irti  ~  0.00069  for  Ar,  Az  =  1/60.  On  further  grid  refinement,  as  illustrated 
by  the  Re  =  2126  case,  ITe\  —  0.00077  for  Ar,  Az  =  1/60  and  Ire\  ~  0.00055  for  Ar, 
Az  -  1/100,  indicating  that  Ar,  Az  =  1/60  constitutes  a  good  compromise  between 
the  level  of  accuracy  and  the  computational  effort  required  to  achieve  it.  The  extent  to 
which  tin-  vorticity  integral  (7)  is  satisfied  is  reduced  as  Re  is  increased  However,  in  the 
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Re  range  considered,  this  trend  is  very  weak,  as  when  Re  =  1492  and  A r,  A z  =  1/60, 
/rel  ~  0.00072,  whereas  for  Re  =  3023  and  Ar,  A z  =  1/60,  ITei  ~  0.00078. 

Although  most  of  the  solutions  presented  in  this  study  reach  a  steady  state  (those 
which  do  not  are  noted),  it  is  important  to  ensure  that  the  model  has  temporal  sis  well  as 
spatial  accurracy.  In  the  parameter  range  where  the  flow  reaches  a  steady  state,  Escudier 
(1984)  reported  that  the  time  taken  to  reach  this  steady  state  from  an  impulsive  start 
was  typically  tens  of  seconds.  In  that  same  parameter  range,  i.e.  1.5  <  H/ R  <  2.0  and 
1000  <  Re  <  1500,  the  numerical  model  reaches  steady  state  for  t  between  500  and  700. 
In  the  case  where  v  =  4.5  x  10 _5m2/s  and  R  =  9.5  x  10_2m,  which  are  the  values  in  the 
experiments  of  Escudier  (1984),  we  find  fi  ~  Re/ 200  and  a  ‘spin-up’  time  to  steady  state 
of  approximately  60  seconds.  The  numerical  result  is  consistent  with  the  experimental 
observations. 

Escudier  found  experimentally  that  tiie  flow  undergoes  a  transition  to  unsteadiness 
at  some  critical  combinations  of  Re  and  H/R  and  he  has  mapped  out  this  transition 
(Figure  1).  At  the  point  where  the  flow  becomes  unsteady,  the  flow  visualization  tech¬ 
nique  is  no  longer  able  to  give  a  clear  picture  of  the  flow  characteristics.  The  present 
numerical  solutions  are  found  to  accurately  predict  those  flows  which  lead  to  a  steady 
state.  This  has  given  us  confidence  to  investigate  the  transition  to  unsteadiness  and  this 
aspect  of  the  flow  is  the  subject  of  a  further  report  which  is  in  preparation. 
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4.  Results  and  Discussion 

Escudier  (1984)  has  mapped  out  a  stability  diagram  in  the  (Re,  H / R)  plane  de¬ 
lineating  the  regions  where  single,  double  and  triple  recirculation  zones  occur,  together 
with  the  boundary  between  steady  and  oscillatory  flows.  This  diagram  is  reproduced  in 
Figure  1.  The  parameter  values  at  which  the  flow  has  been  simulated  in  this  study  are 
also  indicated  in  this  figure.  Escudier  has  made  available  photographs  of  the  recircula¬ 
tion  zones  and  a  number  of  these  are  reproduced  here  for  comparison  with  the  numerical 
solutions.  The  reader  should  be  aware  that  in  the  photographs  of  the  experiment,  radial 
distances  are  reported  to  appear  uniformly  stretched  by  8%  due  to  refraction  at  the 
various  interfaces. 

The  following  is  a  brief  description  of  the  basic  flow  features. 

The  fluid,  which  completely  fills  the  cylindrical  container,  is  initially  at  rest.  At 
t  =  0,  the  bottom  endwall  is  impulsively  set  to  rotate  at  constant  angular  speed  Cl.  An 
Ekman  boundary  layer  develops  on  the  rotating  disk  with  thickness  of  order  (f/fifl2)0  5. 
This  rotating  boundary  layer  then  acts  as  a  centrifugal  pump,  sending  fluid  radially 
outwards  in  a  spiralling  motion  while  ‘sucking’  fluid  into  it  from  above.  This  pumping 
action  of  the  boundary  layer  together  with  the  presence  of  the  cylindrical  wall  at  r  = 
R  sets  up  a  secondary  meridional  circulation.  The  fluid  pumped  out  of  the  Ekman 
layer  spirals  up  the  cylindrical  wall  establishing  a  Stewartson-type  boundary  layer  (see 
Greenspan  1968).  In  time,  fluid  with  angular  momentum  reaches  the  vicinity  of  the 
stationary  top  endwall,  where  it  is  turned  and  advected  towards  the  centre.  It  spirals 
inwards  creating  a  further  boundary  layer  on  the  top  endwall.  On  the  endwall  the  fluid 
separates  at  r  =  0  and  a  concentrated  central  vortex  is  formed  whose  core  size  depends 
on  the  depth  of  the  boundary  layer  from  which  it  emerged.  The  fluid  then  spirals  down 
this  central  vortex  to  be  sucked  back  into  the  Ekman  layer 

The  experimental  studies  of  this  flow  show  that  it  undergoes  a  series  of  bifurcations 
as  the  two  governing  parameters,  Re  and  H/ R,  are  varied.  The  present  study  focusses 
on  the  development  of  the  flow  characteristics  with  Rc  and  H / R  up  to  the  point  where 
the  flow  no  longer  reaches  a  steady  stale.  In  particular,  in  the  region  of  parameter 
space  considered,  i.e.  Re  <  3000  and  // / /?  <  3.5,  the  flow  is  observed  experimentally  to 
remain  axisymmctric  and  laminar,  and  for  the  most  part  to  reach  a  steady  state. 

Figure  1  suggests  that  the  flow  characteristics  depend  critically  on  the  Reynolds 
number,  flic  numerical  solutions  show.  howe\  er,  that  the  basic  dynamics  of  the  central 


vortex  are  inertially  driven  and  that  viscous  shears  are  negligible  outside  recirculation 
zones,  boundary  layers  and  away  from  the  meridional  shear  layer  in  the  bottom  half 
of  the  cylinder  which  is  caused  by  the  strong  turning  at  the  corner  r  =  R,  z  =  0.  A 
comparison  of  the  ‘streamline’  contours,  i.e.  intersections  of  stream  surfaces  with  the 
meridional  plane,  and  contours  of  the  angular  momentum  T  =  rv,  on  the  meridional 
plane  show  that  even  at  the  relatively  low  value  of  Re  —  1000  (see  Figure  2a)  T  is 
constant,  to  high  order,  on  streamlines  outside  boundary  layers.  At  this  value  of  Re ,  the 
curvature  of  the  streamlines  of  the  vortical  fluid  returning  towards  the  lower  rotating 
endwall  does  not  change  sign  and  the  secondary  meridional  flow  consists  of  a  simple 
overturning  motion.  From  the  contour  plots  of  the  azimuthal  vorticity  (Figure  2a(iii)). 
we  find  that  t]  is  positive  in  the  boundary  layer  regions  and  negative  in  the  interior 
flow,  consistent  with  a  deceleration  of  the  meridional  flow  in  the  boundary  layers.  Of 
particular  interest,  one  should  note  the  small  region  of  positive  77  just  above  the  Ekman 
layer  about  the  axis  of  symmetry. 

It  is  interesting  to  examine  the  effect  of  increasing  Re  for  the  case  in  which  ////?  = 
2.5  At  Re  =  1600  a  wavincss  is  evident  in  the  streamtubes  of  the  central  vortex  as 
well  as  in  the  contours  of  T  (see  Figure  2b).  Comparing  the  T  contours  of  the  two  cases 
Re  =  1000  and  Re  =  1600,  it  is  clear  that  a  greater  proportion  of  the  fluid’s  angular 
momentum  has  been  advectcd  towards  r  =  0,  in  the  upper  boundary  layer  region, 
leading  to  the  formation  of  what  appears  to  be  a  weak  centrifugal  wave  in  the  higher 
Re  case.  The  contours  of  tj  at  Re  =  1600  (see  Figure  2b(iii ))  are  qualitatively  similar 
to  those  at  Re  =  1000,  but  now  the  extrema  in  the  values  of  the  azimuthal  vorticity 
have  increased.  The  region  of  positive  azimuthal  vorticity  in  the  interior  of  the  flow  is 
extended  considerably  and  roughly  coincides  with  the  region  of  ‘wavy’  streamtubes. 

The  above  effects  are  further  enhanced  as  Re  is  increased  to  1800.  Here  the  stream- 
tubes  (Figure  2c)  clearly  indicate  the  presence  of  a  centrifugal  wave  with  two  periods 
and  two  bulges  on  the  axis  at  approximately  1/3  and  2/3  of  H / R.  Note  however  that 
the  central  vortex,  ;.s  it  emerges  from  the  endwall  boundary  layer  has  streamtubes  which 
are  almost  parallel  with  a  relatively  large  radial  gradient,  indicating  the  presence  of  a 
relatively  large  axial  velocity  in  this  core  flow. 

An  increase  in  Re  to  1918  crosses  the  boundary  between  no  ‘breakdowns'  and  1 
-  2  ‘break'!'  aviis’  on  Figure  1.  The  amplitude  of  the  waves  has  increased  and  their 
wavelet;  th  A -Teased  (see  Figure  2d).  With  the  increased  amplitud  -  and  decreased 
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wavelength,  the  associated  axial  deceleration  is  large  enough  to  cause  the  flow  to  stagnate 
under  the  crest  of  the  wave.  Within  this  stagnant  region  is  found  a  near  spherical  region 
of  recirculating  fluid  which  is  termed  a  vortex  breakdown  bubble,  i.e.  a  toroidal  vortex 
centered  on  the  axis  of  a  columnar  vortex.  At  this  particular  value  of  H/R  there  are 
two  breakdown  bubbles,  the  downstream  bubble  being  considerably  smaller  and  the 
recirculating  flow  within  it  slower  than  in  the  leading  bubble. 

For  H/R  <  1.95,  both  the  experimental  investigations  of  Escudier  (1984)  and  our 
numerical  calculations  only  identify  single  breakdown  bubbles  and  for  H/R  <  1.2,  no 
breakdown  bubbles  are  found.  Pao  (1970)  observed,  in  a  related  problem  where  the 
cylinder  wall  also  rotates  at  the  same  angular  speed  as  the  disk,  that  for  H/R  =  1.0  and 
Re  >  500  standing  waves  appear  in  the  axial  plane  and  reported  a  similar  dependence  on 
Re  for  these  waves  as  is  found  in  this  study.  His  experimental  investigations  include  cases 
up  to  Re  =  8600  where  the  flow  was  reported  to  be  steady,  axisymmetric  and  laminar. 
However,  he  did  not  find  any  recirculating  breakdown  bubbles  in  his  low  aspect  ratio 
{H/R  =  1.0)  container. 

As  Re  is  increased  to  1942  (Figure  2e),  1996  (Figure  2f)  and  2126  (Figure  2g),  very 
good  agreement  is  found  between  the  calculated  stream  surfaces  and  the  observations  of 
Escudier  (1984).  Some  of  these  experimental  observations  are  included  for  comparison 
with  the  numerical  solutions  in  Figure  3. 

Although  the  length  of  the  recirculation  bubbles  has  grown  larger  with  the  increases 
in  Re,  the  wavelength  of  the  wavy  disturbance  on  the  outer  streamtubes  has  been  re¬ 
duced.  There  remains  a  general  trend  that  as  Re  is  increased,  the  wavelength  is  reduced 
and  the  wave  amplitude  increased.  The  increase  in  the  bubble  diameters  with  Re  is 
indicative  of  the  increased  wave  amplitudes.  The  shortening  of  the  wavelength  is  also 
consistent  with  the  behaviour  of  the  recirculation  bubbles.  As  Re  is  increased,  the  re¬ 
gion  of  near  parallel  flow  just  after  the  vortex  emerges  from  the  boundary  layer  on  the 
stationary  endwall  is  shortened.  This  leads  to  a  migration  of  the  leading  bubble  to¬ 
wards  this  endwall.  Note  also  that  the  distance  between  the  two  bubbles  is  progressively 
reduced  as  Re  is  increased 

The  distribution  of  I'  is  particularly  interesting.  For  fluid  in  the  central  vortex 
that  has  emerged  from  the  upper  boundary  layer,  T  is  conserved  on  the  streamtubes  for 
some  axial  distance  over  which  the  vortex  remains  concentrated  As  the  fluid  spirals 
downstream  past  the  breakdown  region,  there  is  a  change  in  the  angular  momentum. 
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When  the  fluid  is  first  deflected  radially  outwards  around  the  breakdown  bubble,  there  is 
an  increase  in  F  and  when  it  is  converged  again  around  the  downstream  half  of  the  bubble 
it  loses  some  of  this  gain  in  angular  momentum.  The  gradients  in  azimuthal  velocity 
are  large  enough  for  the  viscous  stresses  to  effect  these  changes.  Past  the  breakdown 
regions,  the  angular  momentum  distribution  then  corresponds  to  essentially  solid  body 
rotation  of  the  fluid.  In  essence,  in  these  flows,  the  breakdown  region  is  like  a  transition 
region  from  a  concentrated  vortical  flow  to  solid  body  rotation. 

The  details  of  the  flow  structure  in  the  recirculation  bubble  are  of  particular  interest 
in  this  study,  even  though  their  dynamical  significance  to  the  overall  flow  structure  is 
quite  secondary.  The  main  reason  for  this  interest  is  that  they  provide  an  excellent 
test  for  the  accurracy  of  the  numerical  solutions.  It  is  possible  to  obtain  very  good 
flow  visualizations  of  the  breakdown  bubbles  and  their  structure,  and  correspondingly 
of  critical  points  and  flow  reversals.  These  features  are  much  more  demanding  on  a 
numerical  model  than  the  relatively  simple  structure  of  the  outer  flow. 

Comparing  in  detail  the  structure  of  the  recirculation  regions  as  determined  nu¬ 
merically  and  observed  experimentally  for  Re  =  2126  and  HI R  =  2.5  we  find  very 
good  agreement  between  the  two.  The  structure  of  the  upstream  bubble  is  now  clearly 
interpreted  with  the  aid  of  the  numerical  streamlines.  The  leading  stagnation  point  is 
well  defined  in  the  flow  visualization,  however  the  downstream  stagnation  point  is  not 
as  some  dye  finds  its  way  into  the  bubble  and  some  is  swept  past.  Note  that  the  leading 
dye  streak  is  very  thin  -  this  is  also  reflected  in  the  calculated  streamlines  which  are 
very  close  together  in  the  approach  flow  to  the  stagnation  point.  As  the  fluid  passes 
over  the  bubble  at  maximum  diameter,  the  streamlines  indicate  a  local  acceleration  of 
the  meridional  flow  over  the  ‘obstacle’,  i.e.  bubble.  As  the  fluid  passes  the  bubble,  it 
converges  back  towards  the  axis,  however  the  axial  flow  is  now  considerably  slower  than 
it  was  in  its  approach  to  the  leading  stagnation  point.  This  is  seen  in  the  increased 
distance  between  the  calculated  streamlines  as  well  as  in  the  much  broadened  dyestreak. 
The  dye  photograph  of  the  leading  bubble  shows  a  thin  well  defined  outer  envelop  with 
a  rather  broad,  diffuse  inner  core.  This  is  matched  by  the  computed  streamlines  where 
they  are  contracted  to  the  outer  edge  of  the  bubble  and  spread  out  within  the  interior. 
As  for  the  downstream  bubble,  a  similar  situation  is  found.  This  information  is  not 
evident  from  the  dye  photograph  since  it  appears  that  the  dye  is  residing  principally 
between  a  narrow  band  of  streamlines  which  envelop  the  bubble,  with  virtually  no  dye 
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inside  the  bubble.  The  computed  streamlines  indicate  that  the  recirculating  flow  in  the 
downstream  bubble  is  considerably  slower  than  in  the  upstream  bubble.  Downstream 
of  the  second  bubble,  the  enveloping  dyestreak  is  broad  and  is  tapered  as  it  approaches 
the  Ekman  layer.  This  behaviour  is  mirrored  in  the  computed  streamlines  and  may  be 
interpreted  as  an  axial  acceleration  of  the  flow,  probably  as  a  result  of  the  boundary 
layer  suction. 

This  very  close  agreement  between  experimental  observation  and  numerical  calcula¬ 
tion  of  the  recirculation  zones  has  not  been  achieved  in  every  case,  although  in  all  cases 
the  agreement  has  been  very  good.  The  reason  is  because  the  structure  of  the  recircula¬ 
tion  zones  in  the  vicinity  of  parameter  space  where  a  recirculation  zone  just  emerges  is 
critically  dependent  on  the  value  of  Re  For  example,  the  small  change  in  Re  from  1918 
to  1942  (1.2%)  results  in  a  25%  increase  in  the  diameter  of  the  leading  bubble  and  a  300% 
increase  in  the  diameter  of  the  second  bubble,  whereas  the  outer  flow  remains  virtually 
unchanged.  Hence,  a  small  percentage  error  in  the  estimate  of  Re  can  be  responsible  for 
a  significant  difference  between  the  observed  structure  of  the  recirculation  bubbles  and 
their  computed  structure.  Unfortunately,  no  uncertainty  estimates  on  Re  are  available 
for  the  experimental  observations.  Escudier  reports,  however,  that  the  temperature  of 
the  fluid  was  maintained  at  25°C  ±  0.1"C  and  the  viscosity  of  the  fluid  varied  by  about 
5%  per  °C.  Hence  a  possible  uncertainty  of  0.5%  exists  in  the  estimate  of  Re  due  to  the 
uncertainty  in  viscosity  alone.  For  flows  that  are  not  so  close  to  these  critical  parameter 
regions,  the  structure  of  the  recirculation  bubbles  is  not  so  variable  with  Re  and  a  much 
closer  agreement  between  the  computations  and  the  observations  was  found. 

As  the  Reynolds  number  of  the  flow  increases,  the  downstream  bubble  approaches 
the  tail  end  of  the  upstream  bubble.  Just  beyond  the  value  of  Re  at  which  the  tail 
stagnation  point  of  the  upstream  bubble  coincides  with  the  leading  stagnation  point  of 
the  downstream  bubble,  the  flow  bifurcates.  This  bifurcation  leads  to  a  flow  in  which 
there  are  two  stagnation  points  on  the  axis  of  symmetry  defining  the  head  and  tail  of  the 
recirculation  zone  which  encloses  two  distinct  zones.  These  two  zones  are  connected  by  a 
ring  seen  as  two  saddle-points  either  side  of  the  axis  of  symmetry  in  the  meridional  plane. 
In  the  case  of  Re  =  2494,  H/ R  -  2.5  (see  Figure  2h),  this  flow  pattern  is  steady.  The 
two  critical  streamlines,  i.e.  ip  —  (I  enclosing  the  bubble  and  4’  -  1  609  x  10-6  defining 
the  saddle-points,  are  very  close  together  around  the  outer  envelope  of  the  recirculation 
bubble,  whereas  there  is  a  considerable  gap  between  them  near  the  axis  of  symmetry. 
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This  is  also  reflected  in  the  dye  photograph  (see  Figure  3b)  where  the  dye  depicting  the 
outer  shell  of  the  bubble  is  a  thin  well  defined  shell,  while  it  is  broad  and  diffuse  in  the 
inner  core.  Also,  in  the  vicinity  of  the  saddle-points,  the  dye  lines  are  confused  and  it  is 
difficult  to  resolve  the  flow  from  them. 

When  Re  is  increased  to  2765,  the  computed  flow  is  unsteady  and  its  structure 
oscillates  between  the  structures  typical  of  Re  =  2126  and  Re  =  2494.  Figure  2i  shows 
time-averages  of  r/>,  F  and  tj.  The  development  of  this  unsteadiness  and  its  periodicity  is 
the  subject  of  further  study.  Note,  however,  that  the  outline  of  the  head  of  the  upstream 
recirculation  zone  and  the  outline  of  the  internal  downstream  recirculation  zone  shown 
in  the  ‘snap-shot’  of  the  experimental  flow  (Figure  3c,  found  by  Escudier  (1984)  to  be 
steadily  oscillating)  are  well  represented  by  the  time-averaged  computed  streamlines. 
This  suggests  that  the  dye-lines  of  the  unsteady  flow  retain  some  information  about  the 
time  history  of  the  flow. 

The  flow  develops  into  three  distinct  zones  for  Re  >  2680  when  H/ R  =  3.25  and 
Re  >  3000  for  H/R  =  3.5.  In  all  cases  considered,  the  two  downstream  recirculation 
zones  merge  and  result  in  a  flow  pattern  similar  to  that  described  earlier  for  H/ R  =  2.5 
and  Re  =  2494.  The  upstream  breakdown  region  however  remains  separate  and  its  flow 
structure  is  as  described  for  the  upstream  bubble  when  H/ R  =  2.5  and  Re  =  2126.  The 
numerical  solutions  for  this  triple  breakdown  did  not  reach  steady  state  after  t  =  1000, 
which  corresponds  to  20,000  time  steps.  In  real  time,  this  is  equivalent  to  integrat¬ 
ing  the  system  of  governing  equations,  from  rest,  for  approximately  70  seconds  for  the 
corresponding  experimental  values  of  Escudier  (1984).  The  parameter  values  for  which 
these  solutions  were  obtained  are  very  near  those  delineating  the  unsteady  flow  regions 
according  to  the  stability  diagram  (Figure  1),  where  it  was  observed  experimentally  that 
steady  conditions  are  reached  in  a  very  long  time,  longer  than  10’s  of  seconds.  Perhaps 
if  the  integrations  were  carried  out  for  longer  times  they  would  indeed  reach  steady 
state.  The  oscillations  for  Hj R  =  3.25  are  almost  undetectable,  consisting  of  oscilla¬ 
tions  of  approximately  1%  in  the  values  of  the  streamfunction.  The  oscillations  found 
in  the  solutions  for  H/R  =  3.5  were  considerably  larger,  with  the  middle  breakdown 
zone  oscillating  most.  The  upstream  breakdown  bubble  is  close  to  steady  state,  and 
the  tail  end  of  the  downstream  breakdown  region  is  also  varying  very  slowly.  Figure  4 
are  time-averaged  contours  of  r/>,  T  and  i)  for  Re  =  3061  and  H/R  =  3.5,  which  are  to 
be  compared  with  Figure  5  which  is  a  photograph  of  the  flow  for  the  same  values  of 
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the  parameters.  It  is  interesting  to  note  Escudier’s  (1984)  observation  that  the  central 
recirculation  zones  for  the  H/ R  =  3.5  cases  never  develop  into  well  defined  regions,  and 
this  now  appears  to  be  due  to  these  regions  never  quite  reaching  steady  state. 
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5.  Conclusion 

A  numerical  solution  of  the  axisymmetric  Navier-Stokes  equations  has  been  ob¬ 
tained  and  used  to  examine  the  vortex  breakdown  phenomenon  which  occurs  at  certain 
values  of  the  governing  parameters  for  flows  in  an  enclosed  cylinder  driven  by  a  rotating 
endwall.  The  accuracy  of  the  numerical  solutions  has  been  established  by  varying  the 
grid  resolution  and  satisfying  an  integral  identity  which  exists  for  the  flow.  The  numer¬ 
ical  solutions  are  compared  extensively  with  available  experimental  results,  particularly 
dye-streak  photographs  of  the  flow,  mostly  at  steady  state.  The  extent  of  the  agreement 
between  the  numerical  solutions  and  the  experimental  results  is  critically  examined  and 
is  found  to  be  very  good.  Certain  features  which  were  not  fully  resolved  by  the  ex¬ 
periments  have  now  been  given  plausible  explanations  from  the  numerical  results.  The 
numerical  solutions  also  provide  a  clear  picture  of  the  topology  of  the  flow,  especially  of 
the  structure  of  the  multiple  breakdown  zones. 

In  part  II  of  this  study,  the  numerical  solutions  are  examined  further  and  a  physical 
understanding  of  the  vortex  breakdown  phenomenom  is  gained. 

The  scene  is  also  set  for  future  studies  in  the  unsteady  regime.  Unsteady  flows 
are  generally  less  amenable  to  study  by  flow  visualization  techniques  as  relationships 
between  unsteady  pathlines,  streamlines  and  streaklines  are  not  straightforward.  The 
interpretation  of  streaklines  from  flow  visualization  results  for  example,  requires  more 
care.  For  the  time-dependent  numerical  model  presented  here,  these  distinctions  are 
more  readily  made  and  the  interpretation  of  the  flow  characteristics  more  tractable. 
However,  in  extending  the  present  axisymmetric  model  into  the  unsteady  regime  the 
assumption  of  axial  symmetry  is  in  question  where  there  are  no  experimental  results  to 
support  it.  It  is  not  expected  to  be  valid  as  the  flows  become  increasingly  oscillatory 
and  unstable  as  Re  and  H/R  are  increased. 
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Re 

H/R 

nr 

nz 

i(v) 

/(v+) 

nr) 

1256 

1.5 

40 

60 

-.000386 

.12590 

-.12628 

1256 

1.5 

60 

90 

-.000175 

.12619 

-.12637 

1492 

1.5 

40 

60 

-.000418 

.12888 

-.12929 

1492 

1.5 

60 

90 

-.000188 

.12908 

-.12927 

1854 

1.5 

40 

60 

-.000461 

.13130 

-.13176 

1854 

1.5 

60 

90 

-.000206 

.13417 

-.13438 

1994 

2.5 

60 

150 

-.000237 

.15684 

-.15708 

2126 

2.5 

60 

150 

-.000244 

.15893 

-.15918 

2126 

2.5 

100 

250 

-.000175 

.16041 

-.16058 

2889 

3.5 

60 

210 

-.000277 

.18056 

-.18083 

3023 

3.5 

60 

210 

-.000286 

.18206 

-.18234 

TABLE  1 :  CHARACTERISTICS  OF  THE  VORTICITY  INTEGRAL  (7)  FOR  SELECTED 
CASES. 


FIGURE  1.  STABILITY  BOUNDARIES  FOR  SINGLE,  DOUBLE  AND  TRIPLE 
BREAKDOWNS,  AND  BOUNDARY  BETWEEN  STEADY  AND  OSCILLATORY  FLOW, 
IN  THE  (Re,  H/R)  PLANE  (EMPIRICALLY  DETERMINED  BY  ESCUDIER  (1984)).  THE 
LOCATIONS  IN  PARAMETER  SPACE  WHERE  THE  FLOW  HAS  BEEN  SIMULATED  ARE 
INDICATED  BY  □. 


FIGURE  1 


FIGURE  2:  CONTOURS  OF  (i)  y.  (ii)  T  AND  (Iii)  //IN  THE  MERIDIONAL  PLANE  FOR  H/R  = 

2  5  AND  Re  AS  INDICATED.  THE  CONTOUR  LEVELS  ARE  NON-UNIFORMLY  SPACED, 
WITH  20  POSITIVE  AND  20  NEGATIVE  LEVELS  DETERMINED  BY  CONTOURLEVEL(i) 

=  MA ^variable)  x  (i/20)3  AND  CONTOURLEVEL(i)  =  MIN(i /ariable)  x  (i/20)3 
RESPECTIVELY.  ALL  ARE  PLOTTED  AT  t  =  1000  BY  WHICH  TIME  STEADY  STATE 
FLOW  CONDITIONS  HAVE  BEEN  REACHED,  EXCEPT  FOR  (i)  WHICH  SHOWS  THE  TIME 
AVERAGES  OVER  750  £  t  <  1000  OF  AN  OSCILLATORY  FLOW. 


FIGURE  2a  (i) 


(ii)  n  Re=1000,  H/R  =  2.5 


FIGURE  2a  (iii) 
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FIGURE  2b  (i) 


FIGURE  2b  (ii) 


(iii)  v:  F?e=  1600  H/R  =  2.5 


FIGURE  2b  (in) 


(ii)  r  Re=1800  H/R  =  2.5 


FIGURE  2c  (ii) 


-7.9x10  3 


max.^l.OxlO'6 


FIGURE  2d  (i) 


(ii)  r:  Re=  1918  H/R  =  2.5 


FIGURE  2d  (ii) 


min.  =  -7.8x103 


max.  =  8.8x  =  10'6 


FIGURE  2f  (i) 


(ii)  T:  Re=1994  H/R  =  2.5 


FIGURE  2f  (ii) 


FIGURE  2f  (hi) 


(i)  T:  Re  =  2 126  H/R  =  2.5 


min.=-7.6x10'3  max.  =  3  OxlO-5 


FIGURE  2g  (i) 


(ii)  T:  Re=2126  H/R  =  2.5 


min.=0  max.  =  1.0 


FIGURE  2g  (ii) 


(i)  ‘P:  Re  =  2494  R/E  =  2.5 


min.  =  -7.3x10'3  max.  =  7.4x10-5 


FIGURE  2h  (i) 


(ii)  T :  Re  =  2494  H/R  =  2.5 


t 


i 

< 


FIGURE  2h  (ii) 


min.  =  -4.4  max. =  17. 2 


FIGURE  2h  (in) 


(i)  T:  Re  =  2765  H/R  =  2.5 


min. =-7.0x10-3  max.  =  8.6x10-5 


FIGURE  21  (i) 


FIGURE  3:  VISUALIZATIONS  OF  THE  RECIRCULATION  ZONES  FOR  H/R  =  2.5  AND  Re 
AS  INDICATED.  *  (a)AND(b)  ARE  REPORTED  TO  BE  AT  A  STEADY  STATE,  WHILE  (c)  IS 
A  SNAP-SHOT  OF  A  STEADILY  OSCILLATING  FLOW. 

3(a)  Re  =2126 
3(b)  Re  =  2494 
3(c)  Re  =  2765 


■  A!!  PHOTOGRAPHS  OF  HE  I  XPER1MENTAL  VISUAL!?-  MON  ARE  REPRODU  'ED - 
AIN'  '  PERMISSION  FROM  I.'.  P  I  SCUDIER. 


FIGURE  3c 


FIGURE  4:  CONTOURS  OF  (i)  y,  (ii)  T  AND  (in)  >/  IN  THE  MERIDIONAL  PLANE  FOR  H/R 
3.5  AND  Re  =  3061.  THE  CONTOUR  LEVELS  ARE  NON-UNIFORMLY  SPACED,  WITH 
20  POSITIVE  AND  20  NEGATIVE  LEVELS  DETERMINED  BY  CONTOURLEVEL(i)  = 
MAX(i variable)  x  (i/20)3  AND  CONTOURLEVEL(i)  =  MIN  (variable)  x  (i/20)3 
RESPECTIVELY.  THE  CONTOURS  ARE  OF  TIME  AVERAGES  OVER  750  <  t  <  1000. 


(i)  Re  =  306 1  H/R  =  3.5 


FIGURE  4  (i) 


min. =-4.9 


max.  =  18.2 


FIGURE  4  (Hi) 


FIGURE  5.  VISUALIZATION  OF  THE  RECIRCULATION  ZONES  FOR  H/R  =  3.5  AND 
Re  =  3061 . 


FIGURE  5 
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